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Abstract 

We simulate the dynamics of a quantum dot coupled to the single resonating mode of a metal 
nano-particle. Systems like this are known as metamolecules. In this study, we consider a time- 
dependent driving field acting onto the metamolecule. We use the Heisenberg equations of motion 
for the entire system, while representing the resonating mode in Wigner phase space. A time- 
dependent basis is adopted for the quantum dot. We integrate the dynamics of the metamolecule 
for a range of coupling strengths between the quantum dot and the driving field, while restricting 
the coupling between the quantum dot and the resonant mode to weak values. By monitoring the 
average of the time variation of the energy of the metamolecule model, as well as the coherence 
and the population difference of the quantum dot, we observe distinct non-linear behavior in the 
case of strong coupling to the driving field. 
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I. INTRODUCTION 


Quantum plasmonics is a relatively new area of research that studies the interaction 


of surface plasmons with quantum emitters 
emitter constitute what is known as a metamaterial or metamolecule. Such p 
materials, at variance from those that operate in the microwave regime [6- 


The surface plasmon and the quantum 

asmonic meta- 


8|, require the 


downscaling to nanometers. Metal nano-particles (MNP’s) seem to be an ideal candidate in 


order to build metamaterials at optical frequencies 


Ki. 


They also have the convenient 


property that their resonant frequency can be tuned by changes in their geometries {q]. A 
typical example of such systems is provided by metamolecules comprising MNP’s coupled 
to quantum dots (QD’s) 141. These advances strengthen the need for techniques ca¬ 
pable of simulating the dynamics of MNP’s coupled to QD’s. In the past two decades, 
much progress has been made in the development of methods for the simulation o 
dynamics, involving path integral formulations mean held approximations 


quantum 


classical approximations and surface-hopping schemes 17l422l|. An alternative approach is 


16| . semi- 


provided by the partial Wigner representation of quantum mechanics 


23 


28j. In such a 


representation, exact algorithms for subsystems embedded in harmonic environments can 
be developed (when the coupling to the environment is bilinear) 29 |. 

In this paper, we present a method for simulating, within the partial Wigner representa¬ 
tion of quantum mechanics, a MNP-QD metamolecule subject to an external driving held. 
To this end, we use a piece-wise deterministic algorithm which utilizes a time-dependent 
basis. We study the dynamics of the population diherence of the QD, in the case of weak 
coupling to the MNP, when the QD is subjected to an external driving held of varying 
strengths. By monitoring the average of the time variation of the energy of the meta¬ 
molecule model, as well as the coherence and the population diherence of the quantum dot, 
we observe distinct non-linear behavior in the case of strong coupling to the driving held. For 
validating our approach, we perform a comparison with the results obtained by employing 
a brute-force numerical approach, which uses a discretized phase-space grid and deals with 
the partial diherential equation system associated with the evolution of the density matrix 
in the partial Wigner representation. 

The structure of the paper is as follows. In Sec. [TTl the formalism of quantum dynamics 
in the partial Wigner representation is presented. Section HTTI outlines the generalization of 
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the formalism to a time-dependent Hamiltonian. Section IIVI introduces the model for the 
QD coupled to the MNP. In the same section, a brief outline of the propagation algorithm 
is discussed. In Sec. |Vl the results of the numerical simulations are presented. Finally in 
Sec. rvTl we give our conclusions. Appendix details the adimensional coordinates used in 
our study while App. [B] outlines the phase-space-grid algorithm used to verify the results. 

II. QUANTUM DYNAMICS IN THE PARTIAL WIGNER REPRESENTATION 

Consider a system dehned by the following Hamiltonian operator: 


H = Hs{f, P) + P) +Hc{f,R) + p, R, P, t ), 


( 1 ) 


where S, B and C are subscripts denoting the subsystem, bath and the coupling respectively. 
The Hamiltonian HE{t) describes an external time-dependent held which may interact with 
both the subsystem and bath. The lower case coordinates describe the subsystem degrees 
of freedom, while the upper case coordinates describe the bath. The equation of motion for 
an arbitrary operator x{t) in the Heisenberg picture is written in symplectic form as 


d . i 

= ft 




H 

H m 



t) 



X 


( 2 ) 


where the matrix elements of the symplectic matrix 30(| are dehned as B^- = ej^, with e^j 
being the complete antisymmetric tensor. 

It is assumed that the Hamiltonian of the bath depends on a pair of canonically conjugate 
operators, X = (P, P), and that the coupling Hamiltonian Hq depends only on the position 
coordinates and not on the momenta. The partial Wigner transform for the operator y is 
dehned as 


Xw(-A) 




X 


R + 


( 3 ) 


where X = (P, P) is the phase space point, dehned in terms of the canonically conjugate 
positions and momenta. 

Upon taking the partial Wigner transform of the Heisenberg equation, one obtains the 
Wigner-Heisenberg equation of motion: 


Wt 


x{X,t) 


i 


V 

HwiX) 

h 

Pw(X) Xw{X,t) 



Xw{X, t) 


( 4 ) 
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where 



( 5 ) 


The symbols dk=d /dX^ and dk=d /dX^ denote the operators of derivation with respect 
to the phase-space coordinates acting to the left and right respectively. Snmmation over 
repeated indices is implied. 

The partial Wigner-transformed Hamiltonian takes the form: 


Hw(X, t) = Hs + + HcMR) + t) • 


( 6 ) 


When the Hamiltonians Hb,w, hfc.w and HE,w{X,t) are at most quadratic in R and P, 
the action of the terms in the matrix P are eqnivalent to their linear order Taylor series 
expansion, since any higher order terms acting npon the Hamiltonian will yield zero. Under 
these conditions, algorithms for simnlating exact qnantnm dynamics in the partial Wigner 
representation can be devised. 

The above theory is also applicable in the case when the external field Hamiltonian, 
He, depends npon the bath coordinates, however, in the following we restrict onr stndy to 
sitnations where it depends only upon the subsystem coordinates. 

III. REPRESENTATION IN A TIME-DEPENDENT BASIS 

When considering Eq. (|1]), one can dehne the following time-dependent Hamiltonian 


h^^{R,t) — Hs + + -^E,w(^) ) 


(7) 


where Ub,w(-R) is the potential energy of the bath. A time-dependent basis can be dehned 
in terms of the eigenstates of hw(7?, f): hw(7?, f)|a; i?, f) = Ea{R,t)\a] R,t) . In this basis 
the quantum evolution takes the form 



where to is the initial time, the symbol T denotes time-ordering, and 



( 9 ) 
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In Eq. (j9]), we have introduced the transition operator 


~ + (/3|d )ha/3 , 


( 10 ) 


which arises explicitly from the time dependence of the basis. The operator J^a' 
determines the quantum transitions of the subsystem caused by the interaction with the 


external held. The time-c 
hrst obtained in Refs 


ependent Liouville operator iCaa'piS' is similar in form to those 


^^aa',0IS' f}' T Jaoi',I3jH'i 


( 11 ) 


with = iCjaa'if) + iLaa'{t), and the Bohr frequency given by a;Q,a/(i?, f) = {Ea{R,t) — 
Ea'{R,t))/h. The Liouville operator for the bath degrees of freedom is given by iL^a' = 
{P/M) ■ (d/dR) + (l/2)(F-^(f) + F^) ■ d/dP, where F^{R,t) is a time-dependent Hellman- 
Feynman force for the energy surface Ea{R,t). The quantum transition operator is also 
similar in form to that given in Ref. 32l |: 


Jaa', /3/3'{t) 'Ta^f}{^)^a'f}' T 'T'a'—^l3i{t)SajS , 


( 12 ) 


with 

7 '* 

and AEai 3 {t) = Ea{R, t) — Ei 3 {R, t). In the above, the coupling vector for the time-dependent 
states has been introduced as dafs{R, t) = {a] RR\d /dR\l3] R, t). 

The operator Jaa',/ 3 i 3 '{t) describes the quantum transitions arising from the interaction 
between the system and the bath. If such an interaction is weak, the effect of Jaa',i 3 i 3 '{t) is 
negligible. 


■ dap{R-, t) ( 1 -|- 


p 

M 

— ■ d*^,pi{R, f) ( 1 -h 


1 AEap{t)dap{RR) d 

2 ^■do.^iPR) dP 

lAE^,p>{t)d*^,^,{R,t) d 


^■dlp{RR) 


dP 


(13) 

(14) 


IV. METAMOLECULE MODEL 

The metamolecule model considered in this work comprises a two-level system (the QD) 
coupled to a single resonating mode (RM), and subjected to a time-dependent external 
held. At this point, a single harmonic mode was considered, as the computational resources 
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required by the phase-space grid algorithm, with which results were being compared, rises 
exponentially with the dimension of the bath. In the following, we will use adimensional 
coordinates and parameters; they are expounded in detail in App. 

The QD Hamiltonian is 

Hs = • (15) 

The resonant single-mode Hamiltonian, describing the MNP, is dehned as 

while the coupling to the RM is 

-^c,w = —cRax , (17) 

where c is a coupling constant. The external driving field is represented through the Hamil¬ 
tonian 

^e(^) = 9 cos{uJdt)ax , (18) 


where g denotes the driving strength of the external field, and ojd is the driving frequency. 
The symbols ax and az denote the Pauli matrices. The total Hamiltonian is given by Eq. (| 6 ]). 
The energy eigenvalues of the Hamiltonian hw(-R, t) in Eq. ([7]) are 


Ei^ 2 {R, t) = Vb± \l — +cos"^(ut) 2'yg cos (cat), 


( 19 ) 


where 7 = —cR. In the basis of h^^{R,t), the dynamics of arbitrary quantum operators is 
dehned by Eq. ([S]). One can discretize time and obtain the evolution equation 


Xw'(^) = J^Tjexp 

013' ^ 


Xwito ), 


( 20 ) 


aa' ,13^' 


71 

where Tn = t — R. Using very small time steps and the Dyson identity, one obtains 

X (1 T TnJaa',00' + 'TnJaa',00') ( Xw (O) 


00' n 


( 21 ) 

In the case of weak coupling to the bath, the action of Jaa', 00 ' can be disregarded. For 
Tn = T for every n, one then obtains: 


00' n 


irC^ir) (1 + r|>x^'(to) 


( 22 ) 
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Equation (12^ can be implemented by means of a stochastic algorithm. One can sample with 
probability 1/2 one of the two terms in Jaa',i 3 / 3 ' acting at each time step. For each phase space 
point {R, P), one propagates a single deterministic step, dictated by iC^a' / 3 / 3 '(^)- 
of such a step, the quantum transition, due to the external held, is sampled. Transition 
probabilities can be dehned as: 




^\W)\ 

l + T\{a\/3)\ ■ 


The probability of rejecting the transition will then be given by 


(23) 




1 

l + r|(d|/3)| ■ 


(24) 


For the two-level model that we are studying, the eigenstates can be calculated exactly and 
therefore, so can the transition probabilities in Eqs. fl25]l and fl2T)) . 


V. RESULTS 


The initial state of the system is dehned as 


Pw(-R, P) 



^ X Pb,w(-R, P) , 


(25) 


where the Pb,w(-R, P) is the Wigner function for the bath, given by 

a 


and the matrix on the right hand side of Eq. fl25l) is given in the subsystem basis. The 
average values of an arbitrary phase-space dependent operator, xw{R, PR), are calculated 
as 

{XwiR, P,t)) = Tt' J dRdPp,^{R,P)xw{R-, PR) . (27) 

The partial trace and the evolution of Xw{P, P, t) are calculated in the time-dependent basis 
and with the algorithm sketched in Secs. Hill and HVl 

In order to determine the ehect of the external held upon the metamolecule, we kept 
the values of the system parameters unchanged while varying the coupling strength of the 
external held. The values of the system parameters are /? = 12.5, c = 0.01, U = 0.8, 
oj = 0.5, and ojd = 0.05. These values lead to a weak coupling between the QD and the RM, 



tanh(/3a;/2) 

Pb,w(^,^) =-z-exp 


TT 


2 tanh(/3a;/2) / P^ 

It ^ ^ 


7 









so that the quantum transitions caused by the interaction with the resonant mode can be 
neglected. The values for the coupling strength of the external held used in this study were 
5 ^ = 0.1, 0.3, 0.5, 0.7, 0.9,1.5. However, in order to demonstrate the differences between weak 
and strong coupling, only the results for = 0.1 and g = l.b are shown in the hgures. In the 
case of the piece-wise deterministic algorithm, described in Secs. Hill and HVl a time step of 
r = 0.1 was employed, and a total of 10® trajectories were propagated in each calculation. 
Instead, the phase-space grid algorithm (described in App. requires a smaller time step 
of r = 0.001. The phase-space grid spacing was Ai? = AP = 0.1. 



FIG. 1: Plot of the average coherence, (ux), as a function of time. The solid black circles refer 
to the results obtained with the piece-wise deterministic algorithm while the white triangles refer 
to those of the grid integration. A continuous line connects the data points of the piece-wise 
deterministic algorithm. The values of the system parameters are /? = 12.5, c = 0.01, H = 0.8, 
uj = 0.5, ujd = 0.05 and g = 0.1, corresponding to weak driving field strength. The results of the 
two algorithms are indistinguishable to the human eye. 

In Fig. [H a comparison of the results obtained with the two different algorithms, for the 
average coherence of the quantum dot as a function of time, is shown. With a value of 
5 ^ = 0.1, this calculation corresponds to a weak driving held. The two results agree almost 
exactly, so that they cannot be distinguished by the human eye. Moreover, the error bars 
are negligible. In this weak coupling case, the oscillations remain relatively small, with the 
values of {cTx) ranging between —0.5 and 0.5. A very slow mode of oscillation, with angular 
frequency ~ 0.05 appears to be superimposed to a fast mode with angular frequency ~ 0.81. 
The slow frequency is basically that of the driving held while the slow one corresponds to 
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the tunnel splitting, shifted by a very small amount because of the weak coupling to the 
held. As expected, in the case of weak coupling to the external held it is the tunnel splitting 
that dominates the evolution in time of (cTa,). 
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FIG. 2: Plot of the average coherence, (era,), as a function of time. The solid black circles refer 
to the results obtained with the piece-wise deterministic algorithm while the white triangles refer 
to those of the grid integration. A continuous line connects the data points of the piece-wise 
deterministic algorithm. The values of the system parameters are /3 = 12.5, c = 0.01, If = 0.8, 
a; = 0.5, ujd = 0.05 and g = 1.5, corresponding to strong driving field strength. The results of the 
two algorithms are indistinguishable to the human eye. 

Figure [2] shows the results of the calculations with g = 1.5, corresponding to a strong 
driving held. The results produced by the two algorithms are indistinguishable also in this 
case, with error bars remaining smaller than the points for the entire simulation time. In 
this case, the evolution of the coherence in time displays a non-linear pattern: it starts with 
fast oscillations around (ax) = —0.25 from f = 0 to f ~ 20; within 20 < t < 40 it undergoes 
large oscillations, and then it switches to oscillating fast around the value of (ax) = 0.25; 
it does so until t = 80, when the large oscillations start again. We can therefore conclude 
that the QD switches between two different dynamical regimes, one where the coherence is 
positive, and the other where it is negative. We can associate a period to such a switching 
dynamics, whose numerical value matches that of the driving held. Because of the strong 
coupling, the frequency of the fast oscillations is about 200% greater than that obtained in 
the case of weak coupling to the external held. The inspection of Fig. [2] reveals that the 
dynamics of (ax) is now dominated by the strong coupling to the driving held. 
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FIG. 3: Plot of the average population difference, (cr^j), as a function of time. The solid black 
circles refer to the results obtained with the piece-wise deterministic algorithm while the white 
triangles refer to those of the grid integration. A continuous line connects the data points of the 
piece-wise deterministic algorithm. The values of the system parameters are (3 = 12.5, c = 0.01, 
hi = 0.8, oj = 0.5, Ud = 0.05 and g = 1.5, corresponding to strong driving field strength. The 
results of the two algorithms are indistinguishable to the human eye. 

In Fig. [3l we show the evolution in time of the average of the population difference in 
the strong driving regime. Such a quantity also shows a non-linear pattern, with regions 
of large and fast oscillations separated by regions of slow and small oscillations. At weak 
driving strengths, this behavior is not noticeable, and the oscillations are much smaller. 
This corresponds to a lower percentage of the ensemble of trajectories of the piece-wise 
deterministic algorithm being driven into the excited state by the external held. 

In Fig. m we plot the rate of change of the expectation value of the energy of the QD-RM 
metamolecule in two different cases: g = 0.1 and g = 1.5. For g = 1.5, such a quantity 
shows a non-linear pattern, with regions of large and fast oscillations separated by regions of 
slow and small oscillations, in the same time intervals where {ax) and (a^) do. Instead, for 
g = Q.l the rate of change of the average energy of the metamolecule has signihcantly smaller 
oscillations around the mean and the structured pattern is almost absent. The calculation 
of the rate of change of the energy of RM was also performed, and it was found that the RM 
reacts slowly to the change in energy of the QD. This justihes the neglect of the quantum 
transitions due to the coupling to the RM. 
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FIG. 4: Plot of the rate of change of the expectation value of the Hamiltonian for the metamolecule. 
The main figure displays the results for the strong driving field strength while the inset shows those 
for the weak driving field strength. The continuous line connects the data points. 

VI. CONCLUSIONS 


Employing the partial Wigner representation of quantum mechanics, we have studied the 
dynamics of a model for a quantum dot coupled to a single resonating mode of a metal nano¬ 
particle. We have treated the case in which the Hamiltonian is explicitly time-dependent, 
due to the presence of a driving held directly coupled to the quantum dot. An explicitly 
time-dependent basis has been used for the representation of the equations of motion and 
generalized propagation schemes have been devised both in terms of a piece-wise determinis¬ 
tic algorithm and of a grid-based numerical integration. The results obtained by using these 
two algorithms were compared. We have shown that both schemes of integration produce 
numerically indistinguishable results. However, the piece-wise deterministic algorithm has 
the dehnite advantage of being able to treat systems with a higher number of discrete energy 
levels (such as three- or four-level quantum dots) and many more (hundreds or thousands) 
of resonating modes in an affordable computational time (see, for example. Ref. 29|). 

We have studied the effect of the driving strength of the external held upon the quantum 
dot. By monitoring the average of the time variation of the energy of the metamolecule 
model, as well as the coherence and the population difference of the quantum dot, we 
observe distinct non-linear behavior in the case of strong coupling to the driving held. 

Both the algorithms presented in this work and the results obtained can be considered 
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as a first step toward the development of an effective approach (alternative to the use of 
master equations) for studying plasmonic metamolecules. 


Appendix A: Adimensional coordinates and parameters 


Upon indicating the dimensional coordinates and parameters with a prime, and intro¬ 
ducing the energy scale we have the following dehnitions: 


Q = 

n' 

< ’ 

(Al) 

P = 

p' 

(A2) 


R = 


(A3) 

u = 

OJ' 

< ’ 

(A4) 

c = 

c' 

(A5) 


9 = 

g' 

(A6) 

/S = 

hw'J. 

(A7) 


The symbol M is the inertial parameter of the oscillator, which has been set to unity. Upon 
choosing in such a way that the the frequency u = 0.5 of the oscillator in Eq. ffTOjl 
corresponds to the value u' = 8.9 x 10^^ Hz, which is typical for the dynamics of metal 


nano-particles 


14l |. we obtain a spanned time-scale in our simulations of 5.62 x 10 s and 


an energy variation for the quantum dot, as shown in Fig. 01 of 24.6 meV. 


Appendix B: Phase-space Grid Algorithm 


The equation of the density matrix in the partial Wigner representation is analogous to 
that Given in Eq. dH): 


m 


p{X,t) 


i 


V 

HwiX) 

h 

H^{X) p^{X,t) 



PwiX, t) 


(Bl) 


12 














Using the basis defined by Hs\c() = ^aW) (« = 1, 2), Eq. flBlD can be written as 


dt 


■Pw{R,P,t) 


—tOJaa'Pw ~ PPw ~ \ P^C,wPw ~ Pw -^C,W J 

I {npP^‘' - 

1 ( dKfy^d^' ap^*’'a<w ^ 

2\y dR dP ^ dP dP j ' 


(B2) 


where the frequency Uaa' = {^a — ^a') has been defined. The Liouville operator, L, is 
given by 


j _ p_d _ dHs^w d 

~ dP dP dP ■ 

One can introduce the following frequencies 


(B3) 


= 0 , 0^2 = 0 , 

0^^ = O, = -hvt. (B4) 


The equations of motion for matrix elements of the density operator can be written explicitly 
as: 

—p]^{R,Pp) = ^cR ( 22 lm [p^]) - ^gcos{udt) (2Um [p^]) 

~ -^Pw ~ 2^ ) (1^5) 

P, t) = -iuJ 2 iPw + ^cR (p^ - p^) - cos{uJdt) (pw - Pw) 

- ^Pw ~ (Pw + Pw) . (B6) 

^p^{R,P,t) = -^cR (2Um [p^]) + ^gcos{udt) {2ilm [p^]) 

- ^Pw-|J?(2Re[p^]) . (B7) 

In order to simplify the integration, one can use the definition 

p^“'(X,f) = p^“'(X,f)e-*"““'0 (B8) 
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Hence, using dimensionless coordinates, the equations of motion become 


m 




a 22 
ai’'"' 




-2cR (-Re [q^] sin (a; 2 it) + Im [rj^] cos (a; 2 it)) 

2^1 cos (ujt) (-Re [r]^] sin (wsR) + Im cos (a; 2 R)) 
d 

- c— (Re cos {u 2 it) + Im sin (W 2 R)) , (B9) 

2cR (-Re [q^] sin (£l; 2 R) + Im [q^] cos (W 2 R)) 

251 cos (cut) (-Re [q^] sin (a; 2 R) + Im [q^] cos (a; 2 R)) 

Lq^ - c— (Re [ 77 ^] cos {u 2 it) + Im [q^] sin (W 2 R)) , (BIO) 

-cR (? 7 w - hw) sin(a; 2 R) + fifCOs(a;dt) {q^ - q^) sin(d; 2 R) 

L (Re [q^]) - {q^ + q^) cos{u: 2 it) , (Bll) 

cR - q^) cos{oj 2 it) - gcos{uJdt) {q^ - q^) cos(a) 2 it) 

L (Im [q^]) - (hw + hw) sin(£h 2 R) . (B12) 


Equations fIBOp . fIBlOp . fIBlip and flB12p describe a set of coupled partial differential equa¬ 
tions (PDE’s) which can be solved to obtain the elements of the density matrix as functions 
of time. In order to solve these coupled PDE’s, the numerical integration approach known as 
the method of lines can be employed. The method of lines transforms the PDE’s into a set of 
ordinary differential equations (ODE’s) by using hnite difference approximations for all but 
one of the integration variables. In this case, the phase-space derivatives are approximated, 
while the time variable is left un-approximated. For the types of potentials studied in this 
work, a fourth-order hnite difference approximation for the phase-space derivatives proves 
to be more than sufficient: 


df _ //(X - 2dX) - 8 /(X - dX) + 8 /(X + dX) - /(X + 2dX) 
dX ~ V 12dX 

Once the conversion from PDE’s to ODE’s has been performed, a numerical integration 
method can be utilized. In this work, a Runge-Kutta 5 Cash-Karp method was found to be 
suitable for stable numerical results. In such an approach, the phase-space of the system is 
essentially discretized into a numerical grid, upon which the coupled ODE’s are solved. 
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